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In self-assembly processes, kinetic trapping effects often hinder the formation of thermodynamically stable 
ordered states. In a model of viral capsid assembly and in the phase transformation of a lattice gas, we show 
how simulations in a self-assembling steady state can be used to identify two distinct mechanisms of kinetic 
trapping. We argue that one of these mechanisms can be adequately captured by kinetic rate equations, while 
the other involves a breakdown of theories that rely on cluster size as a reaction coordinate. We discuss how 
these observations might be useful in designing and optimising self-assembly reactions. 



I. INTRODUCTION 

In self-assembly processes, simple components come 
together spontaneously to form ordered products. Such 
processes abound in biology, where the ordered structures 
might be the outer shells of viruses extended one- 
dimensional filaments that make up the cytoskeleton 5 , 
or ordered arrays of proteins on the surface of bacterisP. 
In other areas of nano-science, self-assembled nanostruc- 
tures made from customised DNA oligomers are being 
used to build ever more complex structures 7 , and the 
possibility of tailoring interactions between colloidal par- 
ticles to assemble diverse ordered structures and phases 
is also an area of active researclP. 

This article concentrates on self-assembly processes 
that occur without energy input. From a statistical me- 
chanical viewpoint, such processes may be regarded as 
the relaxation of a system of interacting components to- 
wards their equilibrium state. It is common to draw an 
analogy with phase changes such as crystallisation, as 
discussed long ago by Caspar and Klug in the context of 
viral capsid assembly*! From a theoretical perspective, 
it is natural to separate considerations of self-assembly 
into two parts: Firstly, the thermodynamic (or static) 
problem of determining what equilibrium states can be 
generated by varying the interactions between particles. 
Secondly, there are dynamic questions: how long does it 
take for a system to reach its ordered equilibrium state, 
and how can inter-particle interactions be optimised to 
facilitate rapid and effective assembly? 

This article is concerned with these dynamic questions. 
Even if the equilibrium state of a system is ordered, there 
are many scenarios in which formation of the relevant 
order occurs extremely slowly. In self-assembly, this be- 
havior is often referred to as 'kinetic trapping'; in the 
statistical mechanics of phase change, one more often 
refers to dynamical arrest or to metastable disordered 
states. Practically speaking, these terms indicate that 
self-assembly is rendered ineffective, and that equilibra- 
tion of the system is very slow. 

In recent years, several studies HEDMIU have observed 
that kinetic trapping tends to occur when interparticle 
bonds are very strong, and that assembly is most often ef- 
fective when structures are stabilised by large numbers of 
relatively weak interactions. In such cases, effective self- 
assembly processes are characterised by transient bond 



formation and bond breaking events, leading to the idea 
that microscopic reversibility acts to facilitate effective 
assembly. 

In this article, we use theoretical and computational 
methods from studies of phase change to analyse kinetic 
trapping in self-assembly processes. We concentrate on 
two kinds of kinetic trapping, one of which is familiar 
from classical theories of phase change while the other 
is accompanied by a breakdown of the classical theories, 
and is accompanied by the appearance of many long- 
lived disordered states. We discuss the breakdown of 
these classical theories, and demonstrate that it may be 
identified and characterised through tests of 'local equi- 
libration' assumptions, as proposed by some of We 
illustrate our analysis with a model of viral capsid as- 
sembly and a simple lattice gas model undergoing phase 
separation into dense and dilute phases. 

II. MODELS 

A. Viral capsid model 

We first describe a model for the self-assembly of empty 
icosahedral viral capsids. The model represents capsid 
proteins as rigid bodies ('subunits') with excluded vol- 
ume geometries and orientation-dependent interactions. 
The lowest energy structure is an icosahedral shell con- 
sisting of 20 subunits (details are given in Appendix [A] 
and Fig. [7] as well as in RefP). This model was used 
to examine the assembly of icosahedral viruses around a 
polymer in R efP^ , and is similar to models used by Ra- 
paport et and Nguyen et a/. 15 in simulations of 

empty capsid assembly. Each subunit could correspond 
to a 'capsomer' comprising a trimer of proteins that form 
a T = 1 capsid. 

To simulate the dynamical process of self-assembly, 
we use over-damped Brownian dynamics for the capsid 
subunits as in RefP^, using periodic boundary condi- 
tions and a second order predictor-corrector algorithm^. 
The capsomer subunits have anisotropic translational 
and rotational diffusion constants calculated using Hy- 
drosub7.CP^. To obtain dimensionless units, we rescale 
lengths by <7b , which is the diameter of one of the spheres 
that comprise the excluded volume of the capsomer; 
times are measured in units of to, which is the Brow- 
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FIG. 1. Assembly from a disordered state, (a) Dynamical 
capsid assembly yields in the NVT ensemble. The fraction 
of subunits in well-formed capsids, n ca psid(t) is shown for 
t = 210, OOOto as a function of the binding interaction pa- 
rameter Sb- Snapshots exemplify typical clusters at the cir- 
cled points. Green attractor pseudoatoms are experiencing 
favorable interactions, while gray attractors are not. The size 
of the attractors indicates the length scale of their interac- 
tion. The system contains N — 500 trimer subunits in a box 
of sidelength L = 74<7b (b) Phase change in the lattice gas 
at density p = 0.1. The binodal is located at Sb/T = 1.86. 
The assembly yield is the fraction of particles that have four 
bonds 724 (t). The snapshots show representative configura- 
tions at time t — 10 5 . The dashed line shows the yield that 
would be obtained by equilibrating a very large system: this 
'thermodynamic yield' is monotonic in the bond strength Eb 
while the yield at fixed time is non-monotonic. The lattice 
size is L — 128. 



nian time for a single such sphere. The binding energy 
associated with each interaction site on a subunit is 
and we take Boltzmann's constant /cb = 1 so that the rel- 
evant dimensionless parameter is e^/T. Further details 
of the model are given in Appendix [Al 



B. Viral capsid assembly in the canonical ensemble. 

In Fig. [TJa), we show results from simulations of self- 
assembly at constant particle number, volume, and tem- 
perature (NVT) for various values of the interaction en- 
ergy £b- The initial conditions for the simulations have 
subunits with random positions and orientations. We 
measure the number of perfect capsids in the system 
(a perfect capsid is defined as a cluster with exactly 20 
subunits, each of which has its maximum number of 3 
bonds). We associate the fraction of capsomer subunits 

in perfect capsids n caps ^(t) with the yield of the assembly 

roTi 

process^. 

As anticipated in Section |TJ the yield depends on a 
combination of thermodynamic and kinetic effects. For 
weak bonds (small £b/T), there is little thermodynamic 
drive to assemble and no capsids are formed. For strong 
bonds (large e^/T), the thermodynamic drive to assem- 
bly is strong, but the system is vulnerable to kinetic trap- 
ping and forms disordered clusters of subunits instead of 
perfect capsids. Optimal assembly takes place in an in- 
termediate range e^/T « 4.5 24 . In later sections, we 
will analyse the interplay of thermodynamic and kinetic 
effects in more detail. 



C. Lattice gas model of self-assembly 

We also consider 'self-assembly' in an Ising lattice gas 
containing N particles on a (two-dimensional) square lat- 
tice with V = L 2 sites. Particles may not overlap, so the 
occupancy of site i is £ {0, 1}. Particles on neighbour- 
ing sites form bonds with energy so that the energy of 
a configuration is 
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where the sum runs over (distinct) pairs of nearest neigh- 
bors. Working with a fixed number of particles, the sys- 
tem is unstable to phase separation at low temperatures, 
forming dense (liquid) and dilute (gas) phases. To make 
an analogy with self-assembly, we start with the parti- 
cles in a disordered configuration, and measure the rate 
with which order is formed (see also 14 ). To quantify the 
yield of the assembly process, we measure the number of 
particles that have bonds to all four of their neighbor- 
ing sites: we denote this number by N4 and we write 
77,4 (t) = -^(N^t)) for the fraction of particles with 4 
bonds. 

The model evolves in time according to a Monte Carlo 
(MC) procedure that involves cluster moves, chosen to 
produce trajecto ries that are dynamically realistic, at 
least qualitatively 1 ^^. The method that we use is close 
to that described in RefPn In each MC move, we select 
a seed particle and use it to build a cluster, as follows. 
For each particle bonded to the seed particle, we conduct 
a Monte Carlo trial, adding it to the cluster with prob- 
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This process is then repeated 
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recursively: for those particles that have been added, we 
use the same MC trial to decide whether particles bonded 
to them are added in turn. Taking the resulting cluster, 
we propose a move in a random direction. This move is 
rejected if this proposed move would lead to more than 
one particle on any site. Otherwise, the move is accepted 
with a probability p a = 1/n 2 , where n is the size of the 
cluster to be moved. An MC sweep consists of N moves, 
and time is measured in MC sweeps. The choices of p c 
and p a ensure that the procedure obeys detailed balance 
with respect to a Boltzmann distribution whose energy is 
given by ([!]), and also that large clusters of particles dif- 
fuse freely through the system with a diffusion constant 
consistent with Brownian dynamics D(n) oc 1/n. This 
dynamical scheme represents a schematic description of 
particles with short-ranged attractions moving through 
a solventPU. 

The relevant variables in our NVT simulations of this 
system are the dimensionless bond strength e^/T and 
density p = N/V. The phase behaviour as a function of 
these two parameters is well-known: the system is unsta- 
ble to phase separation at temperatures below the bin- 
odal (that is, when sinh 4 (£: b /2T) > 1/[1 - (2p - l) 8 ]). 
To obtain Fig.[TJb), we initialise particles in random po- 
sitions, and propagate the dynamics. The particles as- 
semble into clusters: for temperatures below the binodal, 
these clusters will grow until their size becomes compara- 
ble to the whole system. However, for the times we con- 
sider, domains are much smaller than the system size, so 
the system is always far from equilibrium. As in the vi- 
ral capsid model, the yield 714 (t) is non-monotonic in the 
bond strength e^/T: the thermodynamic driving force 
to assemble is small when e^/T is small, while kinetic 
trapping sets in for large e\>/T. 



III. STEADY STATE ENSEMBLE - RATE AND 
QUALITY OF ASSEMBLY 

A. Steady state ensemble 

We now discuss a self- assembling steady state, using 
ideas that have been exploited in studies of nucleation 
and phase transformation 28 . Given a self- assembling 
model system such as the capsid or lattice gas model, we 
simulate the dynamics of the system in the usual way, 
except that we periodically remove large clusters of par- 
ticles (subunits) from the system. We refer to clusters 
removed in this way as the products of the self-assembly 
process. The morphologies of the product clusters are 
stored for later analysis, and we then re-introduce free 
particles into the system at random positions so that the 
total number of particles in the system remains constant. 
To make connection with experiment, we imagine a con- 
tinuous assembly process, where free particles (subunits) 
are fed into the system and large assembled products are 
removed, perhaps by exploiting their tendency to sed- 
iment. On starting the system in an initially random 



configuration, it settles down into a time-translationally 
invariant steady state in which product clusters are con- 
tinuously assembling. This feature of the steady state en- 
semble allows time-invariant averages to be taken during 
productive assembly, in contrast to the NVT quenches. 
The criteria for identifying large clusters depend on the 
model being simulated and are described below. 

We use the notation (-) ss for averages within the 
steady state. For example, the average lattice gas en- 
ergy (E(t)) ss is obtained by averaging the energy defined 
in ([TJ at a time t in the steady state regime. We also 
take averages over the product clusters that are formed 
in the steady state. To be precise, after a simulation has 
been in the steady state for a time t b s , let the number 
of product clusters formed in that time be M. Averag- 
ing over many such runs, we obtain the rate of product 
formation, per unit time and unit volume 

R» = -^(M). (2) 

Labelling the product clusters from a given run by an 
index \i = 1,2,...,.M, we may then calculate averages 
over these product clusters, which we denote by (^prod- 
For example, in the lattice gas, if N(fi) is the number 
of particles in cluster \i then (iV(/z)) pro d is obtained by 
averaging this number over all product clusters. 

To make a connection between the steady state en- 
semble and the NVT simulations of Fig. [I] it is useful to 
define the steady state yield 

yss _ gprod x ^ss (3) 

where Q prod is a measure of the 'assembly quality' of the 
products. For example, in the viral capsid model Q prod is 
the fraction of product clusters that are perfect capsids, 
so that the steady state yield Y ss is the production rate 
of perfect capsids. 

B. Viral capsid model in the steady state ensemble 

Results from the steady state ensemble of the viral cap- 
sid model are shown in Fig. [2] In this model, we define 
particles to be connected if they enjoy a bond with in- 
teraction strength U < —bk^T. A cluster of connected 
particles is identified as a product if (i) the cluster has 
either more than 22 particles or it is a perfect capsid, 
and (ii) the cluster has remained bonded with at least 17 
of the same subunits for a time t > 13.3to- We have in 
mind that the products of the assembly process be stable 
long-lived clusters and this second condition limits erro- 
neous identification of weakly-bonded short-lived clusters 
as products. 

In Fig. [2j we show the production rate i? ss , the yield 
F ss , and the product quality Q prod , which is equal to 
the fraction of product clusters that are perfect capsids 
n caps- A s for the NVT simulations of Fig. Ill we observe 
that the yield is nonmonotonic with respect to binding 
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FIG. 2. Steady state ensemble results for the capsid model. We show the yield y ss , the product 'quality' n^ps an d the 
production rate R ss . Four snapshots corresponding to product clusters from the circled perimeter set, 6h = 6.0 are shown on 
top of the plots. The steady state simulations had N = 1000 trimer subunits in a box with sidelength L = 93.23<7b. The red 
arrows indicate the location of hexameric defects discussed in section IV CI 



energy, with an optimum at e^/T « 5. The origin of 
this optimum is a competition between a rate R ss that 
increases on increasing and a quality factor Q(/i)) pro d 
that decreases. (The total production rate increases with 
£b over the whole range considered, although it eventu- 
ally decreases at much higher for reasons discussed in 
section [V] below.) 

In terms of kinetic trapping, we find that for large e^/T 
product clusters are being formed quickly, but these clus- 
ters are of low quality. In later sections we contrast this 
scenario with the 'stalling' or 'starvation' scenarios dis- 
cussed by ZlotnidP^ in the context of viral capsid assem- 
bly. In that case, kinetic trapping appears as a rate R ss 
that decreases sharply as e^jT is large (see also Sec. V| 
below). Fig. [2] shows that this is not the case for the 
steady state ensemble with the parameters we simulate 
here. 



yield nAt) found in Fig. [lj 30 . At p = 0.1 (top panels 
of Fig. [3| the results are similar to those shown for the 
capsid system in Fig. |2j on increasing £h/T, the non- 
monotonic yield arises from a competition between an 
increasing production rate R ss and a decreasing qual- 
ity Q prod . However, in simulations at a lower density 
p = 0.002, the rate R ss is itself non-monotonic. The 
scenario that occurs at low densities is consistent with 
a 'stalling' effedP^, where the system is depleted of free 
particles, leading to slow cluster growth. But it is the rel- 
atively high density (p = 0.1) scenario in the lattice gas 
that mimics the data for the viral capsid model shown in 
Fig. [2] As in that case, kinetic trapping occurs not just 
because of depletion of free particles, but rather from dis- 
ordered large clusters or aggregates, examples of which 
are shown in Fig. [3] 



C. Lattice gas model in the steady state ensemble 

In the lattice gas model, clusters are identified as prod- 
ucts if their size is larger than a maximal cluster size 
n max . We choose n max = 100 although our results depend 
only weakly on n max . For the product quality Q prod , we 
take the fraction of product particles that have 4 bonds, 
calculated as 
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where N^p) is the number of 4-bonded particles in prod- 
uct cluster p. 

Results from the steady state ensemble are shown in 
Fig. [3] for two different densities p. The non-monotonic 
behaviour of the yield Y ss mirrors the behaviour of the 



IV. MEASURES OF CLUSTER EQUILIBRATION 

We now discuss the relation between cluster quality 
Qprod an j a condition that we call 'cluster equilibration'. 
Our idea is that one type of kinetic trapping arises from 
disordered aggregates such as those discussed above, and 
that the importance of these aggregates may be measured 
through deviations from cluster equilibration^. 

For a general definition of cluster equilibration,we char- 
acterise clusters of particles by their size n and by a sec- 
ond index a, /?, 7, . . . that indicates their morphology. If 
Afn,a is the number of clusters with size n and morphol- 
ogy a then our cluster equilibration condition is 



-(E nja — E n ^)/T 



(5) 



where E n ^ a is the energy of a cluster of n particles and 
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FIG. 3. Steady state ensemble in the lattice gas with n m ax = 100. At two different densities, we show the yield y ss , the 
product 'quality' n^ 1 ^ and the 'production rate' R ss . At p = 0.1, we compare with the yield (n±(t)) at a time t — 10 5 after a 
'quench' from a disordered state (data from Fig.[l] rescaled for comparison). In the central panels, we show example 'product' 
clusters, to indicate their morphologies. 



morphology a, and the averages might be taken at a fixed 
time during assembly in the NVT ensemble, or in the 
steady state ensemble. In words, ^ states that: 'for 
clusters of size n, the probabilities of different morpholo- 
gies are Boltzmann-distributed'. It seems natural to in- 
terpret this as a 'cluster equilibration' condition. (If the 
clusters have different excluded volumes one might take 
this into account by replacing the energy with a suitable 
enthalpy, and any internal entropy of the cluster can also 
be incorporated through a cluster free energy.) In the- 
oretical treatments of self-assembly based on rate equa- 
tions or field- theoretic arguments, it is natural to assume 
that ([5| holds (see Sec. |V| below). We now show that de- 
viations from (|5| are signficant throughout the regimes 
where kinetic trapping is important, indicating that such 
deviations must be taken into account in theories of self- 
assembly. 



A. Cluster equilibration in the lattice gas model 

In the self-assembling steady state of the lattice gas 
model, we count the number of four-bonded particles 
within each cluster. We average this quantity over clus- 
ters of a fixed size n, and we denote this average by 
(AT 4 (n)) ss . We emphasise that these are averages over 
clusters in the self-assembling steady state, and not over 
product clusters. It is convenient to compare (iV 4 (7T,)) ss 



with 7Vf s (n), which is the number of four-bonded par- 
ticles in a cluster of size n that minimises the cluster 
energy. We then define 

(An 4 (n)) ss = i(JV 4 (n) - JVf»>ss (6) 

to measure the deviation of the cluster 'quality' from its 
ground state value, normalised by the cluster size n. 

To test the extent of cluster equilibration, we have per- 
formed umbrella sampling, in which we choose a maxi- 
mal cluster size n um b and reject all MC moves that form 
clusters of size bigger than n um b- On propagating the 
dynamics, the system relaxes to a state that satisfies this 
constraint but is otherwise equilibrated, so that we ex- 
pect ^ to holcM 

In the umbrella-sampled ensemble, we again measure 
the number of particles with four bonds and average over 
clusters of size n. The analogue of ([6| within this en- 
semble is (An 4 (n)) umb = £(JV 4 (n) - Nf(n)) umh . Com- 
parison of (An 4 (n)) between ensembles allows a test of 
the cluster equilibration condition: if ([5| holds exactly 
in the self- assembling steady state then (An 4 (n)) ss = 
(An 4 (n)) um b. In Fig. [4j it can be seen that cluster equi- 
libration holds quite accurately at e^/T = 2.5 which 
is close to the maximum of the yield (recall Figs, [ljb) 
and [3]). However, as e^/T increases and assembly qual- 
ity is reduced, a strong departure from cluster equili- 
bration is apparent: we find that (An 4 (n)) ss increases 
while (An 4 (n)) um b decreases. The key point is that 



6 



(a) 

(An 4 (n)> 



0.6 
0.4 
0.2 




- e b =6.0 




=4.0 

b 




- ~ e b =2.5 ^ 




JWvVV\ umb 


steady state 







20 



40 60 

n 



80 



100 



(b) 



M 




10 

cluster size, a? 



FIG. 4. (a) Measurement of cluster equilibration in the lat- 
tice gas model at p = 0.1, by comparison of steady state 
and umbrella-sampled ensembles. The data points show 
(An4(n)) ss while the solid lines show (Ari4(n)) um b. For 
£h/T — 2.5, the system is far from global equilibrium, but 
deviations from the cluster equilibration condition are small 
(compare the black symbols with the solid black lines). As 
£h/T increases, deviations from cluster equilibration increase, 
(b) Similar data for the viral capsid model. The deviation in 
number of bonds between clusters and their ground states is 
shown as a function of the cluster size n for indicated values 
of Eh- All data points correspond to results from the steady 
state ensemble except for the curve with black A symbols, 
which were obtained from umbrella sampling. 



the crossover in and the deviations from 

cluster equilibration occur at similar values of the bond 
strength. Our conclusion is that effective self-assembly 
requires transient bond-breaking processes in order to 
avoid kinetic trapping, and further that these bond- 
breaking processes need to be frequent enough that the 
system is close to the cluster equilibration condition 

Finally, we note that (An^fji)) tends to increase with 
n in a 'sawtooth' fashion. The effect is primarily due 
to the quantity 7Vf s (n) that appears in the definition of 



(An^n)}. As n increases, 7Vf s (n) changes in discrete 
steps of various sizes, depending on the precise nature of 
the cluster ground state. However, there are often a range 
of cluster morphologies with energies close to the ground 
state energy, all of which occur with significant probabil- 
ity in both umbrella-sampled and steady state ensembles. 
The effect of these clusters is that (N^n)} depends more 



smoothly on n than 7Vf s (n), resulting in a sawtooth struc- 
ture in (An^n)). For our purposes, the relevant compar- 
ison is between umbrella-sampled and steady-state data, 
which both exhibit similar n-dependence in this case. 



B. Cluster equilibration in the viral capsid model 

To test cluster equilibration in the viral capsid model, 
we concentrate on the average number of bonds in clus- 
ters of size n, denoted by (B(n)) ss . We compare this 
average with the number of bonds B gs (n) in a cluster 
of size n with minimal energy. In this case the absolute 
deviation from the ground state cluster is of particular 
relevance, since rate equation descriptions of capsid as- 
sembly often consider only the ground state morphology 
for each intermediate size. Therefore we define 



(AB(n)) ss = (B(n) - B gs (n)) s 



(7) 



Note that this deviation is not normalized by the clus- 
ter size n. As for the lattice gas, we perform umbrella 
sampling that prohibits the formation of clusters larger 
than n um b. (Specifically, we use a hybrid Brownian dy- 
namics/Monte Carlo approach where we use a short se- 
quence of unbiased Brownian dynamics steps as a trial 
move, which is rejected if the size of the largest cluster is 
greater than n um b.) 

Results for (AB(n)) are shown in Fig. pj For the 
umbrella-sampled data, we find that (Ai?(n)) um b ~ 
for £b = 4.5: this quantity is similarly small for > 4.5. 
(The apparent deviation from (A£?(n)) um b ~ at n — 18 
in Fig. [4] is likely a result of imperfect equilibration in 
the umbrella sampled simulations). As in the lattice gas 
data, the steady state measurements show that devia- 
tions from cluster equilibration are small near optimal 
assembly, and grow as kinetic trapping sets in and Q prod 
decreases. 

As in the lattice £ as, a sawtooth structure is visible 
in (AB(n)) ss . Here, increasing the cluster size n leads 
to a change of either one or two bonds in B gs (n). As 
(AB(n)) ss deviates from B gs (n), there are several rele- 
vant cluster morphologies in the steady state ensemble 
which average out the step changes that occur in B gs (n): 
typically, the change in (B(n)) ss on increasing n would 
be somewhere between 1 and 2 bonds. The combination 
of discrete changes in B gs (n) and smoother changes in 
(B(n)) ss results in the apparent sawtooth pattern. 

Results in the umbrella sampled ensemble are analysed 
in more detail in Appendix [Cj We find that the free en- 
ergy of a cluster of size n can be obtained by analysing the 
total number of bonds formed together with the entropy 
associated with different ground state morphologies. For 
the purposes of this section, the sawtooth structure in 
Fig. mb) can be attributed to the fact that some cluster 
sizes (n = 5, 8, 10, ... ) have ground states in which every 
capsomer has at least two bonds (see Fig. [5|. For these 
structures B gs (n) is large, but the number of such mor- 
phologies is rather small (see in particular Fig. [8jb)). As 
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FIG. 5. Clusters of sizes (top left to bottom-right) 5, 8, 10, 
and 12 in which every capsomer has at least two bonds. For 
numbers of subunits in between the sizes shown, there are 
no structures in which every capsomer has at least two un- 
strained bonds. This pattern gives rise to the sawtooth form 
for (AB(n)) in Fig. [I] 



cluster equilibration breaks down, the effect on (AB(n)) 
is most pronounced for these cluster sizes, since there are 
many available morphologies with fewer bonds than the 
ground state, and these morphologies tend to form most 
quickly as clusters grow. For other cluster sizes, devia- 
tion from cluster equilibrium are less pronounced, since 
there are diverse ground state morphologies, all of which 
are kinetically accessible. 



V. KINETIC EQUATIONS IN SELF-ASSEMBLY 

In the capsid and lattice gas models, clusters of parti- 
cles grow as assembly takes place. A natural approach is 
therefore to describe this process in terms of kinetic rate 
equations for cluster concentrations. In phase change 
processes, this idea goes back to Becker and Doring 33 , 
and a derivation of this approach from the microscopic 
dynamics of the lattice gas (or Ising) model was consid- 
ered by Binder and Stauffer 3 -^. Similar ideas have been 
developed by Zlotnick and co-workers in order to describe 
viral capsid assembl}^. 

In this section, we show that non-monotonic steady- 
state yields Y ss can be predicted by such equations, but 
we emphasise that these equations fail to capture the 
decreasing quality Q prod that occurs in both capsid and 
lattice gas models. We argue that this failure of kinetic 
rate equations is linked with the breakdown of the cluster 
equilibration condition (J5|. 



A. Equations for cluster growth and self-assembly 

The central idea behind kinetic rate equations is to 
organise configurations of the system according to the 
sizes of the clusters that are present in the system. Let 
N n (t) be the number of clusters of size n, at some time 
t, so that p n (t) = M n (t) /V is the concentration of such 



clusters. For large systems where the various clusters are 
well-mixed and interact through binary collisions, one 
often writes 

^Pn(t) =Y.{ W n-n>,n'Pn-n>(t) - W+ n , Pn (t)] Pn ,(t) 
n' 

+ Yl W n + n',niPn + n>{t) - W^ n \ Pn (t) (8) 
n' 

where the coeffecients W + and W~ are rate constants 
for binary cluster fusion and cluster fission events re- 
spectively. We use a notation where the sums over n' 
are unrestricted, but the coefficients W^ n > are zero for 
n' > n. For a simple description, we may take W+ n , 
and W~ n , to be finite only when n' = 1, recovering the 
classical Becker-Doring equation. 

The restriction to binary collisions may be relaxed 
straightforwardly (see, for eample 28 ) and cases when 
the clusters are not well-mixed can be treated by field- 
theoretic approach^. However, an additional assump- 
tion on writing ^ is that all clusters of size n behave 
statistically identically, regardless of their shape. This 
assumption is tied in with our condition of cluster equi- 



libration above, as discussed in Sec. V C below. 



B. Non-monotonic production rate R in kinetic equations 

The steady state ensemble has a natural realisation in 
terms of these kinetic rate equations. To keep a compact 
notation, we write M = n max — 1 as the size of the largest 
clusters that are not removed as products. We consider 
clusters of sizes n = 1 . . . M, and we restrict ourselves for 
convenience to monomer binding and unbinding. Then, 
for 1 < n < M we have 

d 

— p n (t) = i^pi(t)[pn-iW-PnW]+A n+ ip n+ i(t)-A n p n (t) 

(9) 

For simplicity, we have replaced the n-dependent rate 
constants by a single 'diffusion-limited' rate ZJ^J and A n 
is the rate for unbinding of a monomer from a cluster 
of size n. If the system is allowed to equilibrate, we 
have that p® q = p^ Y\n=2^ Pi 1 / ^n) - In comparing with 
lattice gas or capsid models, we expect monomer binding 
and unbinding rates to be related by detailed balance, as 



D = XrnVe 



Sb/T 



(10) 



where v is an entropic factor associated with bonding, 
with dimensions of volume (specifically, the contribution 
of monomer attractions to the 2nd virial coefficient of the 
system is ve~ £h ^ T ). 

In the assembling steady state, the equations of motion 
for pM (t) and pi (t) are modified to include the removal of 
product clusters: details are given in Appendix [B] The 
total number of particles (subunits) in the system is a 
constant px = ^Z n n Pn- The production rate may also 
be identified as R(t) = Dpi(t)pM(t)- 
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2 3 4 5 6 

e h /T 



FIG. 6. (a) Rate R ss vs ib/T for kinetic rate equations, 
showing non-monotonic behaviour due to 'kinetic trapping' 
in states with many intermediates and few monomers. We 
take M = 50, m* = 10 and the rate is normalised by its value 
as Sb — > oo. (b) The fraction of particles in the assembling 
steady state that are free monomers, further emphasising that 
the small rate for large Sb arises from states with a small num- 
ber of monomers, and hence a small rate of bond formation. 



The simplest case is irreversible binding, where bonds 
never break, so A m = for all m. As shown in Ap- 
pendix [B] the exact result is 



R°° = Dp lPM = 



ADp 2 T 



M 2 (M + 1) 



(11) 



(Within the steady state, we drop all time arguments 
on p n and R.) The signature of kinetic trapping will be 
that introducing some non-zero unbinding rates A n will 
lead to an increase in R (holding pT constant). That 
is, increasing the rate of monomer unbinding increases 
the production rate R. This is the 'stalling' (starvation) 
effect of Zlotnick and co-workers^. 

To observe this effect in the steady state, an essential 
model ingredient is that the unbinding rates A m depend 
on the cluster siz e m. For simplicity and to maintain 
contact with Refs. 29 38 , we suppose that there is a 'critical 
cluster size' m* above which unbinding is slow A m « 
while for small clusters we take a finite value A m = 
A. (The critical cluster size should be interpreted in the 
spirit of classical nucleation theory 3 -^.) 

The production rate R ss depends on m*, n max and 
a dimensionless parameter A/Dpx- This last parameter 
determines the rate of bond-breaking for clusters with 
m < m*: it is convenient to express this as an 'effective 
bond strength' 



e h /T = \og(Dp T /X m ) 



(12) 



Comparison with ([10]) shows that = Sh — Tlog(vpx) 
and is thus a grand free energy with px the concentration 
of a subunit bath. That is, the relevant binding free 
energy depends on the total subunit density as well as 
the bonding parameters and v. The key point is that 
within the rate equation treatment, the full dependence 
of the system on A and pT can be obtained through the 
single parameter Sh/T. (We emphasise however that we 
have assumed that unbinding from large clusters is very 
slow: the rate A in this analysis is the rate of unbinding 
from small clusters.) 



The central result of this analysis is shown in Fig. [6j 
the production rate R shows a non-monotonic depen- 
dence on £b- Since we are working at fixed px, this corre- 
sponds to a non-monotonic dependence on in the cap- 
sid and lattice gas models. Hence, these results qualita- 
tively mirror the behaviour shown in Fig. [3] as well as the 
stalling (or starvation) effects discussed by ZlotniclP^. 
We have verified that the non-monotonicity survives on 
introducing small finite rates for unbinding from large 
clusters (finite A m for m > m*), although a monotonic 
response is recovered if the unbinding rate is completely 
independent of m. 

Physically, the interpretation of this "starvation" 
regime is that a small unbinding rate A acts to reduce 
the concentration of monomers pi since free subunits 
quickly join growing clusters. The production rate is 
R = DpipM so a small concentration pi reduces this 
rate strongly. As the unbinding rate A increases, pi in- 
creases quickly, while the effect on the concentration of 
large clusters pM is much weaker. Thus R increases as A 
is increased, demonstrating that kinetic trapping occurs. 

We note that while we have analysed these kinetic 
equations in the steady state ensemble, similar non- 
monotonic production rates are observed on starting with 
disordered states and waiting for clusters to forrr i 3 * 29 * 38 * 39 !. 



C. Cluster equilibration 

The previous results demonstrate that Eqs. (8| repro- 
duce one feature of the lattice gas and capsid models, 
the nonmonotonic production rate. However, it is clear 
from ([8| that this rate equation approach treats all clus- 
ters of a given size on the same footing. As discussed 
above, these approximations are justified if all clusters of 
a given size behave statistically identically. Classicall}P^, 
the argument supporting this assumption is that large 
clusters are rare, and that transitions between different 
morphologies are rapid compared to collisions between 
clusters. If this separation of time scales holds, one may 
consider each cluster as a separate subsystem, which re- 
laxes quickly to a 'quasiequilibrium' state: the cluster 
equilibration condition ([5| then holds exactly. In prac- 
tice, the condition of cluster equilibration is much weaker 
than the assumption of a clear separation of time scales 
between cluster rearrangement and cluster growth - but 
the results of Sec. HH and [IV] show that it is the clus- 
ter equilibration condition that breaks down as assembly 
quality falls. 

Therefore, when modelling assembly with rate equa- 
tions of this form, there is an implicit assumption that 
cluster equilibration holds, and hence that the assem- 
bly quality Q prod is independent of temperature. From 
Figs. [2] and [3j this assumption is not valid once kinetic 
trapping sets in. Thus, while kinetic rate equations can 
can reproduce a non-monotonic dependence of produc- 
tion rate on bond strength, our results from the steady 
state ensemble show clearly that these equations miss an 
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important part of the story: the decrease of production 
quality as bonds get strong. 

We note that there are two mechanisms by which clus- 
ter equilibration can be violated. In the first, subunits 
form strong interactions with a sub-optimal number of 
partners. In other words, each subunit-subunit interac- 
tion approximately corresponds to a minimum in the in- 
teraction potential, but subunits do not add on to a grow- 
ing cluster in locations that offer the most interaction 
partners. In the second mode of violation, subunits form 
strained bonds which deviate from the ground state of the 
interaction potential. For example, assembling capsids 
frequently form hexameric defects, as illustrated in Fig. [2] 
The first form of cluster equilibration violation can be in- 
corporated into the rate equation approach, at a cost of 
significantly increased computational complexity, if the 
space of all possible cluster configurations can be pre- 
defined, and then relevant cluster configurations can be 
enumerated ahead of time 4 ^ or sampled stochastically 4 ^. 
However, these approaches have not been used to address 
the possibility of defective bonds, for which it is not pos- 
sible to predefine the space of possible cluster configura- 
tions. 



VI. DISCUSSION AND OUTLOOK 

The usefulness of weak interparticle bonds for self- 
assembly lms_been commented on by several au- 
thor p4iiofti4| i9 | 38 | 42 | Thermal fluctuations allow these 
bonds to be broken: we have shown that this effect can 
enhance assembly by increasing the concentration of free 
particles and hence the rate of cluster formation. These 
results are consistent with studies by Zlotnick and co- 
workers. However, our simulations also identify a sec- 
ond mechanism by which weak bonds enhance the as- 
sembly of clusters with a given morphology. Namely, 
bond-breaking processes act to promote cluster equili- 
bration, in the sense of The qualitative importance 
of cluster equilibration was first raised by Whitesides 
and Boncheva^ we have attempted to quantify this idea 
through Eq. 

The importance of kinetic trapping to biological as- 
sembly, and the constraints it places on interactions 
between the constituents, has been vi vidly de mon- 
st rated through exper iment (e.g. Refs.HEHH]) anc [ 
modeling 3 * 4 * 12 * 13 * 15 * 29 * 44 ^. If we are to anticipate the de- 
sign of functionalised particles that self-assemble into or- 
dered structures, the possibility of kinetic trapping must 
surely be taken into account for these systems as well. 
In particular, methods for predicting the "optimal weak- 
ness" of interparticle bonds could streamline the design 
process. Irl^, we proposed that the degree of cluster equi- 
libration (or local equilibration) might be measured us- 
ing fluctuation theorems that couple to the reversibility 
of bond-formation. 

Developments in this direction will be discussed in fu- 
ture publications: here we note that the cluster equilibra- 



tion condition ([5|is weaker than the 'local equilibrium' 
conditions discussed in 4 . For example, ^ may hold even 
in the absence of good- mixing conditions, which lead to 
a deviation from local equilibrium in the sense of* This 
distinction emphasises the point that, while some degrees 
of freedom in out-of-equilibrium systems may be locally 
equilibrated in this sense, other degrees of freedom may 
be far from equilibrium. For example, the recent results 
of Russo and SciortincS^l seem to indicate that density 
fluctuations are much closer to a local equilibrium dis- 
tribution than energy fluctuations. We conjecture that 
the near-local equilibration of density is linked with a 
weak violation of the good-mixing assumption, while the 
energy fluctuations reflect a stronger violation of cluster 
equilibration, in this out-of-equilibrium system. 

More generally, we conclude that our results are en- 
tirely consistent with the general idesP^ that effective 
self-assembly occurs through the reversible formation of 
numerous weak bonds. We believe that statistical me- 
chanical methods can be used to test this idea quanti- 
tatively, with a view to exploiting it in the design and 
control of self-assembly process. In particular, the break- 
down of cluster equilibration when bonds are strong is a 
kinetic effect that is not taken into account in classical 
theories of self-assembly and phase change. We believe 
that the development of quantitative methods for charac- 
terising this effect is a key challenge for theoretical studies 
of self-assembly, and we look forward to further progress 
in this area. 
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Appendix A: Description of the capsid model 

The model subunits are comprised of a set of over- 
lapping spherical 'excluders' that enforce excluded vol- 
ume and spherical 'at tractors' with short-range pair- 
wise, complementary attractions that decorate the bind- 
ing interfaces of the subunit. Each subunit has two 
layers of excluders and attractors. Attractor positions 
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(a) 



(b) 




FIG. 7. The model capsid geometry. (a) Two dimen- 
sional projection of one layer of a model subunit illustrat- 
ing the g eom etry of the capsomer-capsomer pair potential, 
equation ( Al ), with a particular excluder and attractor high- 
lighted from each subunit. The potential is the sum over 
all excluder-excluder and complementary attractor-attractor 
pairs, (b) An example of a well-formed model capsid from a 
simulation trajectory. 



are arranged so that complementary attractors along a 
subunit-subunit interface perfectly overlap in the ground 
state configuration; excluders on either side of the inter- 
face are sepa rated by exactly the cut off of their potential 
(x c , Eq. A2). Subunits have no internal degrees of free- 
dom 



they translate and rotate as rigid bodies. 
The capsid subunits interact through a pairwise poten- 
tial, which can be decomposed into pairwise interactions 
between the elemental building blocks - the excluders 
and attractors. The potential of capsomer subunit z, 
U c &p,h with position Ri, attractor positions {o^} and 
excluder positions {bi} interacting with subunit j is the 
sum of a repulsive potential between every pair of exclud- 
ers and an attractive interaction between complementary 
attractors: 

U cc (Ri, {aj, {bi}, Rj, {bj}, {aj}) = 

^2 A (\Ri + 6? - Rj -b)\, 2 1 /V b , (7b) + 

k,l 
N a 

Xki£hU [\Ri + a\ - Rj - a) \ - 2 1/2 a a , 4<r a , <r £ 

k,l 

(Al) 

where is an adjustable parameter setting the strength 



of the capsomer-capsomer attraction at each attractor 
site, TVb and 7V a are the number of excluders and at- 
tractors respectively, <Jb and <r a are the diameters of the 
excluders and attractors, b^ (a^) is the body-centered 
location of the k th excluder (attractor) on the ith sub- 
unit, xm is 1 if attractors k and / are overlapping in a 
completed capsid (Figure [7]) and otherwise. The di- 
ameter of attractors is set to <r a = for all results 
in this work. 

Lennard- Jones-like potential: 



The function C p is denned as a truncated 



i ((in -(f)- 



■p/2 



• QC q 

: otherwise 

(A2) 

In our dynamical simulations, the capsomer subunits 
have anisotropic translational and rotational diffusion 
constants, calculated as irl^. The unit of time is set by 
the diffusion constant of a single excluder D, and we de- 
fine to = cr^/D. In these dimensionless units, the eigen- 
values of the translational and rotational diffusion ten- 
sors for capsomer subunits are {0.283,0.283,0.197} and 
{0.1906, 0.1906, 0.0984} respectively. 



Appendix B: Production rate within the steady state 

Here we explain how we solved the kinetic equations ([9| 
to obtain the cluster production rate R in the steady 
state ensemble. As discussed in the main text, Eq. ([9| 
with n = M reduces to 

d 

^PmW Dpi(t)[p M -i(t) - pM(t)] - AmPm(^) (Bl) 



and the production rate is 

R(t) = D Pl (t)p M {t). 



(B2) 



For completeness, we also give the equation of motion 
for the monomer concentration p\(t) within the steady 
state, which is 



d_ 

dt 



Pl (t) =MR(t) - 2D Pl (t) 2 + 2\ 2 p 2 (t) 



M 



^2[XnPn(t) ~ Dp 1 (t) Pn _ 1 (t)} (B3) 



In the following, we work in the steady state so we sup- 



press all time depend ence of the p n . Equation (B2) gives 
PM = D Pl /R while ([m]) gives Dp M -i = ^(1 + 
The remaining p n may then be obtained inductively since 
S reduces to D(p n - p n _i) = ^-[A n+ ip n+ i - X n p n ] so 
tnat Pn-i is given in terms of p m with m > n. For 
1 < n < M — 2 we arrive at 



Dp n = 



R 

Pi 



M-l 



7, £ [A 

^ m=n+l 



] (B4) 
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which allows calculation of all of the p n , in terms of R, 
pi and the A n . 

A simple case is when no unbinding takes place, so 
that A m = 0. Then, p n = p\ for all n, and pr = 
M(M+ l)pi/2. Hence the production rate for irreversible 
binding is given by (11). 



(a) 



We now turn to the problem described in the main text, 
where X rn = X for m < m*, with X rn = for m > m*. It 



then follows from (B4) that 



Pn 




, n > m*, 

S(X,m* - n), n < ra*. 



(B5) 



where S(x,n) = (1 — x n+1 )/(l — x) is obtained by sum- 
ming a geometrical progression and A = X/Dpi. We then 
sum over n to obtain px and eliminate R from the result 



using 



R = Dp\/S(k,m* -1) 
[which follows from (|B5|)]. The result is 



(B6) 



(^ r /A)= (M - m * )(M + m * + 1) + /( ^ W * ) (B7) 



with 

/(A,m) = 



2AS(A,ra* - 1) 



rriirn + 1) — 2 m— -\ — 

. ) dX dX 2 



S(A,m) (B8) 



[We used YT=i 

T,?=2 r ( r - l ) xr - 



rx 



n 2 d 2 
dx 2 



S(x,m) and similarly 
S(x,m).] Dimensional analy- 



' dx 



sis shows that the normalised rate R/R°° depends only 
on M, m* and X/Dp^. We therefore fix these parameters 
and solve (B7) numerically for A, obtaining the monomer 
concentration pi = X/(DX). The rate R may then be cal- 
culated from ( |B6[ ) , as shown in Fig. [6] and discussed in 
the main text. 



Appendix C: Binding free energies for the capsid model 

We define the binding free energy (the free energy 
change on adding a capsomer to a cluster of size n) to be 



9h(n) 



-Tin 



P 



n °ss 



Pn-l Pi 



(CI) 



whe re p n is the concentration of clusters of size n (see 
Sec. VA) and c ss is a reference concentration (always 



required when quoting binding free energies). (We take 
&b = 1 throughout this section.) Following RefP^, by 
comparing the size of our capsid to the size of a satellite 



Sa- 



to 



tobacco mosaic virus capsid, we assign c s , 
correspond to 1 M. 

We find that the free energy of dimerization is ap- 
proximately linear in over the range we consider: 




2 3 4 5 6 7 
intermediate size, n 



3 4 5 6 7 
intermediate size, n 



FIG. 8. (a) The binding free energy to add an additional 
subunit is shown as a function of intermediate size for €h — 
4.5. The A symbols denote values computed from umbrella 
sampling simulations, while the □ symbols were calculated 
based on the cluster configurational entropy, as described in 
the text, (b) The change in configurational entropy, As c , 
computed from ground state cluster geometries is shown as a 
function of intermediate size. 



#b (2) ~ -3.5£ b - Ts h + TAs c (2) with the binding en- 
tropy penalty = —10.7 and the configurational en- 
tropy change for dimerization As c (2) = 1.5. Here As c (n) 
is a difference in "configurational entropy", defined as 
As c (n) = ln(f} n /f2 n _i), with ^ n the number of distinct 
ground state cluster configurations with n subunits. (In 
counting distinct configurations, the three edges of each 
capsomer are assumed to be distinguishable, but config- 
urations related by global rotations are not distinct from 
one another. So the number of distinct dimer config- 
urations is 0,2 = (3 2 /2) since there are three possible 
binding sites on each capsomer (hence 3 2 configurations) 
while the factor of 2 accounts for a rotation symmetry 
of the entire dimer.) Note that the value for given in 
RefP^ contains a typo. 

The binding free energy depends on the number of 
contacts that can be formed and the symmetry of the 
ground state complex. To illustrate the latter effect, we 
calculate As c (n) from geometrical considerations. The 
approach follows Zlotnick^ except that we consider all 
possible ground state structures. The resulting configu- 
rational entropy values for intermediates up to size n = 9 
are shown in Fig. [8] (right), where it is evident that clus- 
ters in which every subunit has two or more bonds (e.g. 
n = 5, n = 8) have lower configurational entropy val- 
ues. We discontinued the calculation at n = 10 because 
the number of possible ground state structures becomes 
large, but the trend continues. 

The □ symbols in Fig. [8] are free energies computed 
using the calculated configurational entropy values, with 
the interaction free energy for a single subunit-subunit in- 
terface extracted from gi = gt>(2) + TAs c (2) with gt>(2) 
extracted from the umbrella sampling results, and the 
interaction free energy to form two subunit-subunit in- 
terfaces extracted from g^ — #t>(5) + TAs c (5) . Two 
separate estimates are required since #2 7^ 2#i because 
the binding entropy penalty for forming two bonds is not 
the same as for a single bond. The agreement between 
the extrapolated free energy values and those measured 
from umbrella sampling further illustrates the extent to 



12 



which the system favors ground state configurations at 
equilibrium. 
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